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1 Introduction 

The goal of this paper is to give a mathematical model describing the global behavior of an underground 
waste repository, once the containers start to leak. The purpose of such a global model is to be used 
for the full held simulations used in safety assessements. The physical situation can be described as an 
array made of high number of leaking modules inside a thin low permeable layer (e.g. clay), included 
between two bigger layers with higher permeability (e.g. limestone or marl). The pollutant is trans- 
ported both by the convection produced by the water flowing slowly (creeping flow)through the rocks 
and by the diffusion coming from the dilution in the water. The leaking last the all period of time ]0, t m [, 
that is small compared to the millions of years over which convection and diffusion are active. In a real 
repository there is a pressure drop producing the flow crossing a large number of disposal modules where 
each module includes several containers. Herein, for simplicity, the repository consists of a set of modules 
lying on a hypersurface E and we represent the leaking of a disposal module by a localized density source 
inside the domain or by a hole in the domain with a given flux on its boundary. Moreover, without lost 
of generality we assume the convection velocity field to be given . According to the test case ||, the 
typical size of a module is a hundred of meters for the width, a kilometer for the length and five meters 
for the height. The distance between two modules is also of order 100 meters and the low permeable layer 
(the clay layer), in which the repository is embedded, has respectively a height and a length of order 150 
and 3000 meters. Since there is a large number of modules, each of them with a small size compared 
to the layers size, see figure 1, a direct numerical simulations of the full field, based on a microscopic 
model taking in account all the detail, is unrealistic. Considering the ratio between the width of a single 
module I and the layer length L, which is of order 1/30, as a small parameter, e, in the microscopic 
model, then the modules, have a height of order e 2 , and are now inmbedded in a layer of thickness e. 
The study of the renormalized model behavior, as e tends to , by means of the homogenization method 
and boundary layers, gives an asymptotic model which could be used as a repository global model for 
numerical simulations. 
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Figure 1: Three layers of soil containing alveoli 

We use methods similar to those applied to the fluid flow through a sieve in §, § or §. Similar 
stationary problem with zero source term (i.e. $ = 0) was treated in 



2 Setting the problem 

2.1 The Geometry 

Let ft c R™ be a bounded domain. Let £ C R n_1 x {0} be such that £ c ft. Let A c R" _1 be a 
periodic set obtained by periodic repetition of bounded closed set M CC] - 1/2, 1/2 [ n_1 . More precisely 

A= (J M a , 

qEZ"" 1 

where M a = a + M. For small parameter e -C 1 and (3 > 1 we define B e = (e A fl S) X ] ~ e",e" [ 
(in situation described in the introduction = 2). We denote by J(e) = {a € Z n_1 ; e M a / U } 
and by T E a = d(e M a x]- e f \ e f3 [) . Finally ft e = ft \ S e , ftj = (ft \ S e ) x ]0, T[ , ft T = ft x ]0, T[ and 

r e = 9s e , rj = r e x]o,T[. 

2.2 The Equations 

Let $ G L°°([0, T]) be the function describing the time behaviour of an alveolus. As mentioned before 
it has a compact support [0,t m ] c]0, T[. Let A = > 0, with r being the half life of the radioactive 
element, and let ipo € H (ft E ) the initial concentration of the radioactive material in the soil (typically 
equal to zero). The diffusion is described by A £ L°°(R; R" x ") a positive definite matrix function. Since 
layers of soil involved in our model have different properties, we assume that 



MVn) = 




for \y n \ < h 
for \y n \ > h 



Now we write our diffusion matrix in the form A £ (x n ) = A(^). In the above described situation the low 
permeable layer has a hight of 150 m meaning that, in that case, h = 3/2. We have the same situation 
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with the convection velocity v £ C([0, T]; iJ 1 (R x 0) n ). The dependence on t/„ is similar as in the case 
of diffusion matrix: 

V 1 (x,t) , for \y n \ < h 
v 2 (x,i) , for \y n \> h 



v(x,y n ,t) = 



For simplicity, we assume that the last component v n does not depend on y n . Next, we suppose that 
div x v = 0, in order to have the divergence free convection velocity. Finally we pose v e (x, t) — v(x, — , t). 
At last, we define the porosity of the medium as 



w(t/„) 

and we put 



lo 1 , for \y n \ < h 
w 2 , for \y n \ > h 



u> £ (x n ) = u>(x n /e) . 

The process is governed by the following convection-diffusion type equation: 

u £ ^ - div(A £ Vip £ ) + (v e • V)ip e + \lu £ tp £ = in (1) 

9? e (0, x) = ipa(x) ie!! E (2) 

n • (A £ V<p £ - v e ip e ) = $(*) on (3) 

We also need to impose some boundary condition on the exterior boundary S = dQ. Let S = <Si U 1S2, 
where Si are disjoint and connected parts of S. We impose 

<p £ = on Si (4) 

n • A £ V^ e - (v £ • n) (p e = on S* 2 . (5) 

3 A priori estimates 

The main result of this section is: 

Proposition 1 Let (p e be a unique solution of ^)-^)- Then there exists a constant C > independent 
of e such that 

keU~ ( n ? ) < C (6) 

|Ve|l,2(0,T;Hi(f2e)) < C ( 7 ) 



Proof. The estimate @ is the consequence of the maximum principle. To prove (Q) we use <^ e as the 
test function in (Q)-@- We obtain 



L 2 (n c 



(A £ V(p e |V^e)L2 (f2 T ) +A I 



<P.\h 
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< 

L 2 (n e ) 



<C(l + |^ e | ffW) ) . □ 
4 Weak convergence 

Our solution <p e is defined on variable domain fij. To use the weak convergence methods, we extend it 
to whole domain f2 T preserving the estimates @ , (0) . In the sequel we assume that ip e is extended using 
the results from Q and we denote that extension by the same symbol. Due to the proposition ^, we can 
conclude that there exists some ip G L 2 (0, T; H 1 ^)) n L°°(f2 T ) such that (up to a subsequence) 

tp E -±p weak* in L°°(Q T ) (8) 

Vtp £ Vv? weakly in L 2 (ft T ) . (9) 

The main goal of this section is to identify the limit ip . We prove that: 

Theorem 1 The limit function ip is the unique solution of the problem 

uj 2 |p - div{A 2 Vp) + (v 2 • V)(p + Xuj 2 p = in h T = (Q\£) x ]0,T[ (10) 

^(x,o) = <^o(x) x g n = o\s (ii) 

= on Sx (12) 
n • A 2 \7p - (v 2 • n) <p = on S 2 (13) 
M = , [e„ • A 2 Vp - (v 2 • e n ) tp]=2§ \M\ on E , (14) 

where [w](x') = w(x', 0+) — iu(:c', 0— ) , x 1 = (x\, . . . , x n _\) denotes the jump over E and \M\ denotes the 
area of M . 

Proof. Let ip G C°°([0, T];Cg°(n)) be such that ip{ ■ ,T) = . Using ^ as the test function in 
we get 



= - 



+ ( tu £ \p e ip+ ( $ Y [ ^ 



r-r 



Passage to the limit for the first four integrals is straightforward. For the last integral we have 
i>(x, t) dT\ = ty{xlt) + 0(e) ) |rf| = ^(xf , t) 2 \M\ e^ 1 + 0{e n+ ^ 2 ) , 

where x\ = ( {xf)', 0) is an arbitrary point from e M a x {0}. But then 

J2 [ i>{x,t)dT\^2\M\ f ip(x',0,t)dx' , 
ieJ(e) Jv l Js 
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where x' = (aci, . . . , € R™ 1 . □ 



Remark 1 In fact we did not use the periodicity of distribution of alveoli. The same proof holds in case 
if each alveolus is randomly placed in a mash of an e net. The alveoli do not even need to have the same 
shape, only the areas of their surfaces need to be equal. 

5 Asymptotic expansion 

The above weak limit describes the global long time behaviour of the process in case when the flux $ is 
not too large . However if we need more accurate information on the behaviour in the near field (i.e. in 
vicinity of S), more precise asymptotics is needed. 

To avoid cumbersome computations we simplify the geometry by assuming that fi =] — L/2, L/2[ n . We 
denote then 

E =}- L/2,L/2[ n - 1 x{0} 
S+ =}- L/2,L/2[ n - 1 x{L/2} 
S~ =} -L/2,L/2["- 1 x{-L/2} . 

We also suppose that the alveoli are rectangular (which is true in real- world situation). More precisely 
we take 

ra-1 

M = JJ ] - TO,, TO, [ , 1/2 > TO, > . 

t=l 

We impose the Dirichlet condition on the bottom and Neumann condition on the top of the domain: 



(p £ = on S~ (15) 
n • A £ Vip E - (v e • n) ip e = on 5+ . (16) 

On the lateral boundary we impose the periodicity condition in order to avoid the lateral boundary layer. 
More precisely, we impose 

ip £ is L — periodic in x = (x\, . . . , x n -i) (17) 

We also assume that L/e G N. For the sake of compatibility, we suppose that given data w 1 ,^ 2 and <p 
are L-periodic in x' . 

We expect some fast changes of solution in vicinity of containers. Therefore, in that region, we in- 
troduce the fast variable y = x/e to describe that behaviour. Far from the sources we expect <p e , the 
solution of -(pi), ( p6| ) and (17), to behave almost like our weak limit tp, which, in this case, satisfies 
(p0|), ( p"l| ) , (|l4| ) plus the conditions 

ip = on r (18) 
n • A 2 V<p - (v 2 • n) ip = on S + . (19) 
tp is L — periodic in x' = {x\, . . . , x n _ i) . (20) 
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That suggests the use of method of matched asymptotic expansions (see e.g. g). 
We separate the domain in three parts separated by S+ , Ej : 

(1+ =} -L/2,L/2[ n - 1 x]delog(l/e),L/2[ 
Q- =} - L/2, L/2 r 1 x ] - L/2, -de log(l/e) [ 
G e =] - L/2, L/2 ["^x] -delog(l/ e ),<felog(l/e)[ . 

Constant d > is to be determined later in order to minimize the error of approximation. As suggested, 
in ilf, we approximate ip e by 99° that satisfies the equation (0) and boundary conditions (|l2|), (|l^) as 
well as the initial condition (|ll|). In G e we look for the asymptotic expansion of ip e , in the form 



QtrP 

p e (x,t) » ^(x,t) + e [ X£ fe ( a; /e) ^(z, t) + u> e (a;/e) $(t) ] + 

+s 2 [ x f(x/s) ^^ + w»{x/e) dip l^ t] v(x,t), + $(*) z k e {x/e) v(x,t) k ] + ■■■ (21) 

Here and in the sequel we assume the summation from 1 to n over the repeating index. The function 
tp® mimics the behaviour of ip but has two close jumps in stead of one. In fact, that suggests that more 
accurate approximation of the real situation would be to have two jumps of the flux; one just above and 
another just below the array of alveoli. However taking the weak limit smears those two jumps into one. 
Namely <p® is defined by 

dia° 

u e -j£ - div (A e Vp° £ ) + (v £ • V )<p° + A uj e p° £ = in flj = (fl\(S+ U E~) ) x ]0, T[ 

tp° e (x,o) = cp (x) xen e = n\(E+ u sj) 

<p° = on 5+ (22) 
n • A 2 V<p° - (v 2 -n)<p° = on 5" 

[^] = , [e„ • A 2 V^ - (v • e„) <p°] = - i * |0P e | on E± , 



is i — periodic in a/ = (2:1, . . . , x n -x) , 



with 



E± =] - L/2, L/2 [x{±de Iog(l/e)} , 

[to] (a:') = w(x',de log(l / e) +)- w(x' ,ds bg(l/e) -) a;' G E+ 
M(a:') = w(x', -de log(l/e) +) - w(x' , -de log(l/e) -) x' G E~ 

The functions , X^, and w e arc the solutions of the auxiliary problems of the stationary diffusion 
type posed in an infinite strip 

& = 1/2,1/2^ xR)\P e , 

with 

P e =Mx]-eP-\s i3 - 1 \ . 



G 



First two problems read 



Figure 2: Strip Q E 



-div(AVXe) = in G e 

n • AV( X k £ +y k ) = on dP e (23) 
Xe is 1 - periodic in y' = (y x , . . . , y n _ x ) 
lim Vx k e = ■ 



-div(A Vw £ ) = in £ £ 

n • A Vw £ = 1 on dP e (24) 
w e is 1 - periodic in y' = (t/i, . . . , y n _i) 

lim AVw £ (y) = |9Pe| e„ . (25) 

Solvability of problem ( g^ ) is classical (see e.g. [|| or Js) ). Due to the symmetry of the domain Q e 
we obviously have that 

x k M = x k A-y) ■ (26) 

Furthermore, there exists a constant c k (e), such that 

\x k £ -c k (s)\ m(lynl>s) <Ce-™ , (27) 

for some C, t > . 

Since is determined up to a constant we may assume in the sequel that c k (e) = 0. 

Remark 2 In general we should have two stabilisation constants c^.(e) at ±oo. Since ( p6| ) holds those 
two constants are equal. In case of general P e , considered in the first part, this seams not to be the case. 

The problem ( p4[ ) does not admit a solution with decaying gradient (due to the source term on dP e ). 
Therefore we have imposed (|2^). To see its behaviour at oo we need to cut-off that boundary condition 
first. To do so we take a cut-off function 

for - 1/2 < y n < 1/2 
C(ffn) = { 1 for \Vn\ > 1 

smooth otherwise 



Now we take 

tt( it \ = —r(i, \ ( 4. 
The function v e (y) — w £ (y) — n(y n ) satisfies the problem 
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-div(AVu e ) = (A nn 7r')' in G s 

n • A Vw e = on dP £ (28) 
v £ is 1 - periodic in = (y 1} . . . , y n -i) 
lim Vi> e = 

y n — too 

and it is obviously pair v s (y) — v £ (—y). 

Since the right-hand side is compactly supported, reasoning as in the case of (23), we conclude that such 
problem admits a unique (up to a constant) solution satisfying 

\v £ ~c(s)\ H i {lVnl>s) <Ce- rs , (29) 

where the constant c(e) can be chosen equal to zero. Therefore the asymptotic behaviour of w e , for large 
\Vn\ is 

We{y) ~ -(A^„) _1 |y„| - |<9P E | + exponentially decaying part . (30) 
The auxiliary problems for the second corrector are as follows 

-div (A V X f") = A tt ^ + ^-(A U X?) in Q e 

n • AVxf = on dP £ (31) 
X e e rn is 1 - periodic in y' = (y 1: . . . ,y n -i) 
lim Vxf l = . 

-div(AV<) = M in g £ 

oyt 

n ■ A Vw\ 3 = on dP £ (32) 
w % J is 1 - periodic in y' = (yi, ■ ■ ■ ,y n -i) 
lim Vwj J = . 

-div(AVz £ fc ) = in G £ 

dy k 

n • AVz e fc = on dP e (33) 
is 1 - periodic in y' = (y 1; . . . , 
lim Vz^ = , k^n 

lim (AWz n + ^(Al n )- 1 \dP e \ \y n \)=Q . 
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As before we conclude that Xs an d , k ^ n can be chosen to decay exponentially towards 

zero as y n — > ±00, while 

z e(v) ~ -7 \ dp e\ ( A D~ 2 \Vn\ Vn + exponentially decaying part . 



Now we still have the term ew{x/e) «3> in the inner approximation that hasn't been matched by our 
exterior approximation. Taking into account ( |30| ) we need to patch our outer approximation with a term 
of the following form 

d e log(l/e) <pl . 

At the same time we will correct the flux jump created by z™ . To do that we define the second corrector 



¥>l(x,o) = o x e o e 

^ = on S + 

n • A 2 V^ - (v 2 ■ n) y\ = on S" 

fo 1 ] = T \ $ (AD" 1 |flP e | , [e„ • (A 2 V^J - v 2 ^)] = (AD" 1 |9P £ | on £± , 

(fil is L — periodic in a;' = (x%, . . . , a; n -i) ; 
For expansion (|^) we can prove the following error estimate 

Theorem 2 Let m = /or n = 3 and m < 2 /or n = 2. Lef d > 2. There exists a constant C > 
independent on e, such that 

\<p e -F E \ L2{0tT . HHBe)) <C(elog(l/e)r ■ (34) 

w/iere 

{¥>°(x,t) +de log(l/e) <^(x,t) , m Q± 
^(x, t)+rfe Iog(l/e) ^(x, t) + s[ X k e (x/e) + d log(l/e) e ^ )(x, i) + u> £ (x/e) $(t) ] + 

B £ = fi\(E+US-) 

withY.f =] — L/2, L/2[ n ~ 1 x {±delog(l/e)}. Furthermore the same estimate holds in L°°(0, T; L 2 (f2 £ )) 
norm. 

Proof. We divide the domain in three parts, ttf, f2~, G £ . 
Let i? £ = ip e — F E . In we have 



J ' div(A 2 Vi? £ ) + (v 2 • V)R S + Xlo 2 i? £ = 



dt 



In G £ the function R £ satisfies 

. dR F 



div(A £ Vi? e ) + (v £ • V)R £ + Xto e R e = E £ , 



at 

[R £ ] = 0(e 2 log 2 (l/e)) , [(A 2 WR £ - i? £ v 2 ) • e„ ] = 0{e d ) 



with \E e \ Lx ,r Ge \ < C e log(l/e) . Furthermore, on T,f- we have jumps 



Now the result follows by the standard a priori estimate. □ 

It should be noticed that the terms in our expansion F e still depend on e implicitly. However it is clear 
that: 

Lemma 1 



\<P° - <p\l'(o,T;H1(Q)) < Cy/de log(l/e) 

IVeli2(o,riHHn±) - C^fde \og{\/e) . 

The proof is straightforward and follows by deducing two problems and estimating the remainder. 
We have the following consequences of theorem 0: 



Corollary 1 Let 



He(x) = \ + e[ X k E {x/e) + w £ (x/e) *(f) ] , in Gj 



Then 

Furthermore 



\<p e - H £ \ L 2 {0iT . h i {Be)} < C(elog(l/e) ) 



3/2 



\<Pe ~ <p\l*(o,T;HH(i s )) < C y/ e log(l/e) 
\<Pe ~ ^U»(o,r ; L»(n.)) < C(e log(l/e)) 3 / 2 

For the auxiliary problems we have: 
Lemma 2 

|VCtf-x*)M«w<Cfe (/, - 1)/a 

|V(xr fe -X mfc )lL 2(ee) KCeW-V' 2 

/2 



\V(w £ ~w)\ L 2 m <Ce^/ 2 
\W(wT k -w mk )\ L 2 {gB) <Ce^/ 2 , 
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where \ iX i w i z i w are the solutions of corresponding auxiliary problems posed on 



g =] - 1/2, l/2[ X R\(Af X {0}) 



. M x {0} 



Figure 2: Strip Q 



Remark 3 It should be noticed that the domain Q is not locally placed on one side of its boundary. 
Furthermore, all the problems ( pj[ ) defining \ k f or k n have only trivial solutions. On the contrary \ n 
admits a nontrivial solution and the corresponding auxiliary problem now reads 



-div(AVx n ) = in Q 
d(x n + Vn) 



d%jk 



on M x {±0} 



(35) 



X n is I- periodic in y' = {y x , . . . , y n -\) 



lim Vx" = 



— div(A Vra) =0 in Q 

TA„ fc = 1 on M x {±0} (36) 

oyk 

w is 1 — periodic in y' = (j/i, . . . , y n -i) 
lim AVw = T\M\ e n . 

y n — >-zbco 

6 Conclusion 

The expansion (|2l| ) clearly points out two important terms in the asymptotic behaviour of ip £ , zero order 
term ip° e and first order term e w £ (x/e) In a real- life situation, that we are trying to model, the 
containers are leaking intensively for very short time. During that time $ is large and the second order 
term e w e (x/e) $ dominates the behaviour of the solution ip e (despite of e multiplying it). Indeed, the 
typical diffusion coefficient in a low permeable layer (clay) is small, compared to the one in the rest of 
the domain (limestone). Thus, at the begining of the process, diffusion arround sources is slow. After a 
short period of time $ vanishes and it is only then that the diffusion becomes dominant, i.e. ip® becomes 
the most important term. 
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